home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / fermi_dirac.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  35.1 KB  |  1,634 lines
  1. /* specfunc/fermi_dirac.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_exp.h>
  26. #include <gsl/gsl_sf_gamma.h>
  27. #include <gsl/gsl_sf_hyperg.h>
  28. #include <gsl/gsl_sf_pow_int.h>
  29. #include <gsl/gsl_sf_zeta.h>
  30. #include <gsl/gsl_sf_fermi_dirac.h>
  31.  
  32. #include "error.h"
  33.  
  34. #include "chebyshev.h"
  35. #include "cheb_eval.c"
  36.  
  37. #define locEPS  (1000.0*GSL_DBL_EPSILON)
  38.  
  39.  
  40. /* Chebyshev fit for F_{1}(t);  -1 < t < 1, -1 < x < 1
  41.  */
  42. static double fd_1_a_data[22] = {
  43.   1.8949340668482264365,
  44.   0.7237719066890052793,
  45.   0.1250000000000000000,
  46.   0.0101065196435973942,
  47.   0.0,
  48.  -0.0000600615242174119,
  49.   0.0,
  50.   6.816528764623e-7,
  51.   0.0,
  52.  -9.5895779195e-9,
  53.   0.0,
  54.   1.515104135e-10,
  55.   0.0,
  56.  -2.5785616e-12,
  57.   0.0,
  58.   4.62270e-14,
  59.   0.0,
  60.  -8.612e-16,
  61.   0.0,
  62.   1.65e-17,
  63.   0.0,
  64.  -3.e-19
  65. };
  66. static cheb_series fd_1_a_cs = {
  67.   fd_1_a_data,
  68.   21,
  69.   -1, 1,
  70.   12
  71. };
  72.  
  73.  
  74. /* Chebyshev fit for F_{1}(3/2(t+1) + 1);  -1 < t < 1, 1 < x < 4
  75.  */
  76. static double fd_1_b_data[22] = {
  77.   10.409136795234611872,
  78.   3.899445098225161947,
  79.   0.513510935510521222,
  80.   0.010618736770218426,
  81.  -0.001584468020659694,
  82.   0.000146139297161640,
  83.  -1.408095734499e-6,
  84.  -2.177993899484e-6,
  85.   3.91423660640e-7,
  86.  -2.3860262660e-8,
  87.  -4.138309573e-9,
  88.   1.283965236e-9,
  89.  -1.39695990e-10,
  90.  -4.907743e-12,
  91.   4.399878e-12,
  92.  -7.17291e-13,
  93.   2.4320e-14,
  94.   1.4230e-14,
  95.  -3.446e-15,
  96.   2.93e-16,
  97.   3.7e-17,
  98.  -1.6e-17
  99. };
  100. static cheb_series fd_1_b_cs = {
  101.   fd_1_b_data,
  102.   21,
  103.   -1, 1,
  104.   11
  105. };
  106.  
  107.  
  108. /* Chebyshev fit for F_{1}(3(t+1) + 4);  -1 < t < 1, 4 < x < 10
  109.  */
  110. static double fd_1_c_data[23] = {
  111.   56.78099449124299762,
  112.   21.00718468237668011,
  113.   2.24592457063193457,
  114.   0.00173793640425994,
  115.  -0.00058716468739423,
  116.   0.00016306958492437,
  117.  -0.00003817425583020,
  118.   7.64527252009e-6,
  119.  -1.31348500162e-6,
  120.   1.9000646056e-7,
  121.  -2.141328223e-8,
  122.   1.23906372e-9,
  123.   2.1848049e-10,
  124.  -1.0134282e-10,
  125.   2.484728e-11,
  126.  -4.73067e-12,
  127.   7.3555e-13,
  128.  -8.740e-14,
  129.   4.85e-15,
  130.   1.23e-15,
  131.  -5.6e-16,
  132.   1.4e-16,
  133.  -3.e-17
  134. };
  135. static cheb_series fd_1_c_cs = {
  136.   fd_1_c_data,
  137.   22,
  138.   -1, 1,
  139.   13
  140. };
  141.  
  142.  
  143. /* Chebyshev fit for F_{1}(x) / x^2
  144.  * 10 < x < 30 
  145.  * -1 < t < 1
  146.  * t = 1/10 (x-10) - 1 = x/10 - 2
  147.  * x = 10(t+2)
  148.  */
  149. static double fd_1_d_data[30] = {
  150.   1.0126626021151374442,
  151.  -0.0063312525536433793,
  152.   0.0024837319237084326,
  153.  -0.0008764333697726109,
  154.   0.0002913344438921266,
  155.  -0.0000931877907705692,
  156.   0.0000290151342040275,
  157.  -8.8548707259955e-6,
  158.   2.6603474114517e-6,
  159.  -7.891415690452e-7,
  160.   2.315730237195e-7,
  161.  -6.73179452963e-8,
  162.   1.94048035606e-8,
  163.  -5.5507129189e-9,
  164.   1.5766090896e-9,
  165.  -4.449310875e-10,
  166.   1.248292745e-10,
  167.  -3.48392894e-11,
  168.   9.6791550e-12,
  169.  -2.6786240e-12,
  170.   7.388852e-13,
  171.  -2.032828e-13,
  172.   5.58115e-14,
  173.  -1.52987e-14,
  174.   4.1886e-15,
  175.  -1.1458e-15,
  176.   3.132e-16,
  177.  -8.56e-17,
  178.   2.33e-17,
  179.  -5.9e-18
  180. };
  181. static cheb_series fd_1_d_cs = {
  182.   fd_1_d_data,
  183.   29,
  184.   -1, 1,
  185.   14
  186. };
  187.  
  188.  
  189. /* Chebyshev fit for F_{1}(x) / x^2
  190.  * 30 < x < Inf
  191.  * -1 < t < 1
  192.  * t = 60/x - 1
  193.  * x = 60/(t+1)
  194.  */
  195. static double fd_1_e_data[10] = {
  196.   1.0013707783890401683,
  197.   0.0009138522593601060,
  198.   0.0002284630648400133,
  199.  -1.57e-17,
  200.  -1.27e-17,
  201.  -9.7e-18,
  202.  -6.9e-18,
  203.  -4.6e-18,
  204.  -2.9e-18,
  205.  -1.7e-18
  206. };
  207. static cheb_series fd_1_e_cs = {
  208.   fd_1_e_data,
  209.   9,
  210.   -1, 1,
  211.   4
  212. };
  213.  
  214.  
  215. /* Chebyshev fit for F_{2}(t);  -1 < t < 1, -1 < x < 1
  216.  */
  217. static double fd_2_a_data[21] = {
  218.   2.1573661917148458336,
  219.   0.8849670334241132182,
  220.   0.1784163467613519713,
  221.   0.0208333333333333333,
  222.   0.0012708226459768508,
  223.   0.0,
  224.  -5.0619314244895e-6,
  225.   0.0,
  226.   4.32026533989e-8,
  227.   0.0,
  228.  -4.870544166e-10,
  229.   0.0,
  230.   6.4203740e-12,
  231.   0.0,
  232.  -9.37424e-14,
  233.   0.0,
  234.   1.4715e-15,
  235.   0.0,
  236.  -2.44e-17,
  237.   0.0,
  238.   4.e-19
  239. };
  240. static cheb_series fd_2_a_cs = {
  241.   fd_2_a_data,
  242.   20,
  243.   -1, 1,
  244.   12
  245. };
  246.  
  247.  
  248. /* Chebyshev fit for F_{2}(3/2(t+1) + 1);  -1 < t < 1, 1 < x < 4
  249.  */
  250. static double fd_2_b_data[22] = {
  251.   16.508258811798623599,
  252.   7.421719394793067988,
  253.   1.458309885545603821,
  254.   0.128773850882795229,
  255.   0.001963612026198147,
  256.  -0.000237458988738779,
  257.   0.000018539661382641,
  258.  -1.92805649479e-7,
  259.  -2.01950028452e-7,
  260.   3.2963497518e-8,
  261.  -1.885817092e-9,
  262.  -2.72632744e-10,
  263.   8.0554561e-11,
  264.  -8.313223e-12,
  265.  -2.24489e-13,
  266.   2.18778e-13,
  267.  -3.4290e-14,
  268.   1.225e-15,
  269.   5.81e-16,
  270.  -1.37e-16,
  271.   1.2e-17,
  272.   1.e-18
  273. };
  274. static cheb_series fd_2_b_cs = {
  275.   fd_2_b_data,
  276.   21,
  277.   -1, 1,
  278.   12
  279. };
  280.  
  281.  
  282. /* Chebyshev fit for F_{1}(3(t+1) + 4);  -1 < t < 1, 4 < x < 10
  283.  */
  284. static double fd_2_c_data[20] = {
  285.   168.87129776686440711,
  286.   81.80260488091659458,
  287.   15.75408505947931513,
  288.   1.12325586765966440,
  289.   0.00059057505725084,
  290.  -0.00016469712946921,
  291.   0.00003885607810107,
  292.  -7.89873660613e-6,
  293.   1.39786238616e-6,
  294.  -2.1534528656e-7,
  295.   2.831510953e-8,
  296.  -2.94978583e-9,
  297.   1.6755082e-10,
  298.   2.234229e-11,
  299.  -1.035130e-11,
  300.   2.41117e-12,
  301.  -4.3531e-13,
  302.   6.447e-14,
  303.  -7.39e-15,
  304.   4.3e-16
  305. };
  306. static cheb_series fd_2_c_cs = {
  307.   fd_2_c_data,
  308.   19,
  309.   -1, 1,
  310.   12
  311. };
  312.  
  313.  
  314. /* Chebyshev fit for F_{1}(x) / x^3
  315.  * 10 < x < 30 
  316.  * -1 < t < 1
  317.  * t = 1/10 (x-10) - 1 = x/10 - 2
  318.  * x = 10(t+2)
  319.  */
  320. static double fd_2_d_data[30] = {
  321.   0.3459960518965277589,
  322.  -0.00633136397691958024,
  323.   0.00248382959047594408,
  324.  -0.00087651191884005114,
  325.   0.00029139255351719932,
  326.  -0.00009322746111846199,
  327.   0.00002904021914564786,
  328.  -8.86962264810663e-6,
  329.   2.66844972574613e-6,
  330.  -7.9331564996004e-7,
  331.   2.3359868615516e-7,
  332.  -6.824790880436e-8,
  333.   1.981036528154e-8,
  334.  -5.71940426300e-9,
  335.   1.64379426579e-9,
  336.  -4.7064937566e-10,
  337.   1.3432614122e-10,
  338.  -3.823400534e-11,
  339.   1.085771994e-11,
  340.  -3.07727465e-12,
  341.   8.7064848e-13,
  342.  -2.4595431e-13,
  343.   6.938531e-14,
  344.  -1.954939e-14,
  345.   5.50162e-15,
  346.  -1.54657e-15,
  347.   4.3429e-16,
  348.  -1.2178e-16,
  349.   3.394e-17,
  350.  -8.81e-18
  351. };
  352. static cheb_series fd_2_d_cs = {
  353.   fd_2_d_data,
  354.   29,
  355.   -1, 1,
  356.   14
  357. };
  358.  
  359.  
  360. /* Chebyshev fit for F_{2}(x) / x^3
  361.  * 30 < x < Inf
  362.  * -1 < t < 1
  363.  * t = 60/x - 1
  364.  * x = 60/(t+1)
  365.  */
  366. static double fd_2_e_data[4] = {
  367.   0.3347041117223735227,
  368.   0.00091385225936012645,
  369.   0.00022846306484003205,
  370.   5.2e-19
  371. };
  372. static cheb_series fd_2_e_cs = {
  373.   fd_2_e_data,
  374.   3,
  375.   -1, 1,
  376.   3
  377. };
  378.  
  379.  
  380. /* Chebyshev fit for F_{-1/2}(t);  -1 < t < 1, -1 < x < 1
  381.  */
  382. static double fd_mhalf_a_data[20] = {
  383.   1.2663290042859741974,
  384.   0.3697876251911153071,
  385.   0.0278131011214405055,
  386.  -0.0033332848565672007,
  387.  -0.0004438108265412038,
  388.   0.0000616495177243839,
  389.   8.7589611449897e-6,
  390.  -1.2622936986172e-6,
  391.  -1.837464037221e-7,
  392.   2.69495091400e-8,
  393.   3.9760866257e-9,
  394.  -5.894468795e-10,
  395.  -8.77321638e-11,
  396.   1.31016571e-11,
  397.   1.9621619e-12,
  398.  -2.945887e-13,
  399.  -4.43234e-14,
  400.   6.6816e-15,
  401.   1.0084e-15,
  402.  -1.561e-16
  403. };
  404. static cheb_series fd_mhalf_a_cs = {
  405.   fd_mhalf_a_data,
  406.   19,
  407.   -1, 1,
  408.   12
  409. };
  410.  
  411.  
  412. /* Chebyshev fit for F_{-1/2}(3/2(t+1) + 1);  -1 < t < 1, 1 < x < 4
  413.  */
  414. static double fd_mhalf_b_data[20] = {
  415.   3.270796131942071484,
  416.   0.5809004935853417887,
  417.  -0.0299313438794694987,
  418.  -0.0013287935412612198,
  419.   0.0009910221228704198,
  420.  -0.0001690954939688554,
  421.   6.5955849946915e-6,
  422.   3.5953966033618e-6,
  423.  -9.430672023181e-7,
  424.   8.75773958291e-8,
  425.   1.06247652607e-8,
  426.  -4.9587006215e-9,
  427.   7.160432795e-10,
  428.   4.5072219e-12,
  429.  -2.3695425e-11,
  430.   4.9122208e-12,
  431.  -2.905277e-13,
  432.  -9.59291e-14,
  433.   3.00028e-14,
  434.  -3.4970e-15
  435. };
  436. static cheb_series fd_mhalf_b_cs = {
  437.   fd_mhalf_b_data,
  438.   19,
  439.   -1, 1,
  440.   12
  441. };
  442.  
  443.  
  444. /* Chebyshev fit for F_{-1/2}(3(t+1) + 4);  -1 < t < 1, 4 < x < 10
  445.  */
  446. static double fd_mhalf_c_data[25] = {
  447.   5.828283273430595507,
  448.   0.677521118293264655,
  449.  -0.043946248736481554,
  450.   0.005825595781828244,
  451.  -0.000864858907380668,
  452.   0.000110017890076539,
  453.  -6.973305225404e-6,
  454.  -1.716267414672e-6,
  455.   8.59811582041e-7,
  456.  -2.33066786976e-7,
  457.   4.8503191159e-8,
  458.  -8.130620247e-9,
  459.   1.021068250e-9,
  460.  -5.3188423e-11,
  461.  -1.9430559e-11,
  462.   8.750506e-12,
  463.  -2.324897e-12,
  464.   4.83102e-13,
  465.  -8.1207e-14,
  466.   1.0132e-14,
  467.  -4.64e-16,
  468.  -2.24e-16,
  469.   9.7e-17,
  470.  -2.6e-17,
  471.   5.e-18
  472. };
  473. static cheb_series fd_mhalf_c_cs = {
  474.   fd_mhalf_c_data,
  475.   24,
  476.   -1, 1,
  477.   13
  478. };
  479.  
  480.  
  481. /* Chebyshev fit for F_{-1/2}(x) / x^(1/2)
  482.  * 10 < x < 30 
  483.  * -1 < t < 1
  484.  * t = 1/10 (x-10) - 1 = x/10 - 2
  485.  */
  486. static double fd_mhalf_d_data[30] = {
  487.   2.2530744202862438709,
  488.   0.0018745152720114692,
  489.  -0.0007550198497498903,
  490.   0.0002759818676644382,
  491.  -0.0000959406283465913,
  492.   0.0000324056855537065,
  493.  -0.0000107462396145761,
  494.   3.5126865219224e-6,
  495.  -1.1313072730092e-6,
  496.   3.577454162766e-7,
  497.  -1.104926666238e-7,
  498.   3.31304165692e-8,
  499.  -9.5837381008e-9,
  500.   2.6575790141e-9,
  501.  -7.015201447e-10,
  502.   1.747111336e-10,
  503.  -4.04909605e-11,
  504.   8.5104999e-12,
  505.  -1.5261885e-12,
  506.   1.876851e-13,
  507.   1.00574e-14,
  508.  -1.82002e-14,
  509.   8.6634e-15,
  510.  -3.2058e-15,
  511.   1.0572e-15,
  512.  -3.259e-16,
  513.   9.60e-17,
  514.  -2.74e-17,
  515.   7.6e-18,
  516.  -1.9e-18
  517. };
  518. static cheb_series fd_mhalf_d_cs = {
  519.   fd_mhalf_d_data,
  520.   29,
  521.   -1, 1,
  522.   15
  523. };
  524.  
  525.  
  526. /* Chebyshev fit for F_{1/2}(t);  -1 < t < 1, -1 < x < 1
  527.  */
  528. static double fd_half_a_data[23] = {
  529.   1.7177138871306189157,
  530.   0.6192579515822668460,
  531.   0.0932802275119206269,
  532.   0.0047094853246636182,
  533.  -0.0004243667967864481,
  534.  -0.0000452569787686193,
  535.   5.2426509519168e-6,
  536.   6.387648249080e-7,
  537.  -8.05777004848e-8,
  538.  -1.04290272415e-8,
  539.   1.3769478010e-9,
  540.   1.847190359e-10,
  541.  -2.51061890e-11,
  542.  -3.4497818e-12,
  543.   4.784373e-13,
  544.   6.68828e-14,
  545.  -9.4147e-15,
  546.  -1.3333e-15,
  547.   1.898e-16,
  548.   2.72e-17,
  549.  -3.9e-18,
  550.  -6.e-19,
  551.   1.e-19
  552. };
  553. static cheb_series fd_half_a_cs = {
  554.   fd_half_a_data,
  555.   22,
  556.   -1, 1,
  557.   11
  558. };
  559.  
  560.  
  561. /* Chebyshev fit for F_{1/2}(3/2(t+1) + 1);  -1 < t < 1, 1 < x < 4
  562.  */
  563. static double fd_half_b_data[20] = {
  564.   7.651013792074984027,
  565.   2.475545606866155737,
  566.   0.218335982672476128,
  567.  -0.007730591500584980,
  568.  -0.000217443383867318,
  569.   0.000147663980681359,
  570.  -0.000021586361321527,
  571.   8.07712735394e-7,
  572.   3.28858050706e-7,
  573.  -7.9474330632e-8,
  574.   6.940207234e-9,
  575.   6.75594681e-10,
  576.  -3.10200490e-10,
  577.   4.2677233e-11,
  578.  -2.1696e-14,
  579.  -1.170245e-12,
  580.   2.34757e-13,
  581.  -1.4139e-14,
  582.  -3.864e-15,
  583.   1.202e-15
  584. };
  585. static cheb_series fd_half_b_cs = {
  586.   fd_half_b_data,
  587.   19,
  588.   -1, 1,
  589.   12
  590. };
  591.  
  592.  
  593. /* Chebyshev fit for F_{1/2}(3(t+1) + 4);  -1 < t < 1, 4 < x < 10
  594.  */
  595. static double fd_half_c_data[23] = {
  596.   29.584339348839816528,
  597.   8.808344283250615592,
  598.   0.503771641883577308,
  599.  -0.021540694914550443,
  600.   0.002143341709406890,
  601.  -0.000257365680646579,
  602.   0.000027933539372803,
  603.  -1.678525030167e-6,
  604.  -2.78100117693e-7,
  605.   1.35218065147e-7,
  606.  -3.3740425009e-8,
  607.   6.474834942e-9,
  608.  -1.009678978e-9,
  609.   1.20057555e-10,
  610.  -6.636314e-12,
  611.  -1.710566e-12,
  612.   7.75069e-13,
  613.  -1.97973e-13,
  614.   3.9414e-14,
  615.  -6.374e-15,
  616.   7.77e-16,
  617.  -4.0e-17,
  618.  -1.4e-17
  619. };
  620. static cheb_series fd_half_c_cs = {
  621.   fd_half_c_data,
  622.   22,
  623.   -1, 1,
  624.   13
  625. };
  626.  
  627.  
  628. /* Chebyshev fit for F_{1/2}(x) / x^(3/2)
  629.  * 10 < x < 30 
  630.  * -1 < t < 1
  631.  * t = 1/10 (x-10) - 1 = x/10 - 2
  632.  */
  633. static double fd_half_d_data[30] = {
  634.   1.5116909434145508537,
  635.  -0.0036043405371630468,
  636.   0.0014207743256393359,
  637.  -0.0005045399052400260,
  638.   0.0001690758006957347,
  639.  -0.0000546305872688307,
  640.   0.0000172223228484571,
  641.  -5.3352603788706e-6,
  642.   1.6315287543662e-6,
  643.  -4.939021084898e-7,
  644.   1.482515450316e-7,
  645.  -4.41552276226e-8,
  646.   1.30503160961e-8,
  647.  -3.8262599802e-9,
  648.   1.1123226976e-9,
  649.  -3.204765534e-10,
  650.   9.14870489e-11,
  651.  -2.58778946e-11,
  652.   7.2550731e-12,
  653.  -2.0172226e-12,
  654.   5.566891e-13,
  655.  -1.526247e-13,
  656.   4.16121e-14,
  657.  -1.12933e-14,
  658.   3.0537e-15,
  659.  -8.234e-16,
  660.   2.215e-16,
  661.  -5.95e-17,
  662.   1.59e-17,
  663.  -4.0e-18
  664. };
  665. static cheb_series fd_half_d_cs = {
  666.   fd_half_d_data,
  667.   29,
  668.   -1, 1,
  669.   15
  670. };
  671.  
  672.  
  673.  
  674. /* Chebyshev fit for F_{3/2}(t);  -1 < t < 1, -1 < x < 1
  675.  */
  676. static double fd_3half_a_data[20] = {
  677.   2.0404775940601704976,
  678.   0.8122168298093491444,
  679.   0.1536371165644008069,
  680.   0.0156174323847845125,
  681.   0.0005943427879290297,
  682.  -0.0000429609447738365,
  683.  -3.8246452994606e-6,
  684.   3.802306180287e-7,
  685.   4.05746157593e-8,
  686.  -4.5530360159e-9,
  687.  -5.306873139e-10,
  688.   6.37297268e-11,
  689.   7.8403674e-12,
  690.  -9.840241e-13,
  691.  -1.255952e-13,
  692.   1.62617e-14,
  693.   2.1318e-15,
  694.  -2.825e-16,
  695.  -3.78e-17,
  696.   5.1e-18
  697. };
  698. static cheb_series fd_3half_a_cs = {
  699.   fd_3half_a_data,
  700.   19,
  701.   -1, 1,
  702.   11
  703. };
  704.  
  705.  
  706. /* Chebyshev fit for F_{3/2}(3/2(t+1) + 1);  -1 < t < 1, 1 < x < 4
  707.  */
  708. static double fd_3half_b_data[22] = {
  709.   13.403206654624176674,
  710.   5.574508357051880924,
  711.   0.931228574387527769,
  712.   0.054638356514085862,
  713.  -0.001477172902737439,
  714.  -0.000029378553381869,
  715.   0.000018357033493246,
  716.  -2.348059218454e-6,
  717.   8.3173787440e-8,
  718.   2.6826486956e-8,
  719.  -6.011244398e-9,
  720.   4.94345981e-10,
  721.   3.9557340e-11,
  722.  -1.7894930e-11,
  723.   2.348972e-12,
  724.  -1.2823e-14,
  725.  -5.4192e-14,
  726.   1.0527e-14,
  727.  -6.39e-16,
  728.  -1.47e-16,
  729.   4.5e-17,
  730.  -5.e-18
  731. };
  732. static cheb_series fd_3half_b_cs = {
  733.   fd_3half_b_data,
  734.   21,
  735.   -1, 1,
  736.   12
  737. };
  738.  
  739.  
  740. /* Chebyshev fit for F_{3/2}(3(t+1) + 4);  -1 < t < 1, 4 < x < 10
  741.  */
  742. static double fd_3half_c_data[21] = {
  743.   101.03685253378877642,
  744.   43.62085156043435883,
  745.   6.62241373362387453,
  746.   0.25081415008708521,
  747.  -0.00798124846271395,
  748.   0.00063462245101023,
  749.  -0.00006392178890410,
  750.   6.04535131939e-6,
  751.  -3.4007683037e-7,
  752.  -4.072661545e-8,
  753.   1.931148453e-8,
  754.  -4.46328355e-9,
  755.   7.9434717e-10,
  756.  -1.1573569e-10,
  757.   1.304658e-11,
  758.  -7.4114e-13,
  759.  -1.4181e-13,
  760.   6.491e-14,
  761.  -1.597e-14,
  762.   3.05e-15,
  763.  -4.8e-16
  764. };
  765. static cheb_series fd_3half_c_cs = {
  766.   fd_3half_c_data,
  767.   20,
  768.   -1, 1,
  769.   12
  770. };
  771.  
  772.  
  773. /* Chebyshev fit for F_{3/2}(x) / x^(5/2)
  774.  * 10 < x < 30 
  775.  * -1 < t < 1
  776.  * t = 1/10 (x-10) - 1 = x/10 - 2
  777.  */
  778. static double fd_3half_d_data[25] = {
  779.   0.6160645215171852381,
  780.  -0.0071239478492671463,
  781.   0.0027906866139659846,
  782.  -0.0009829521424317718,
  783.   0.0003260229808519545,
  784.  -0.0001040160912910890,
  785.   0.0000322931223232439,
  786.  -9.8243506588102e-6,
  787.   2.9420132351277e-6,
  788.  -8.699154670418e-7,
  789.   2.545460071999e-7,
  790.  -7.38305056331e-8,
  791.   2.12545670310e-8,
  792.  -6.0796532462e-9,
  793.   1.7294556741e-9,
  794.  -4.896540687e-10,
  795.   1.380786037e-10,
  796.  -3.88057305e-11,
  797.   1.08753212e-11,
  798.  -3.0407308e-12,
  799.   8.485626e-13,
  800.  -2.364275e-13,
  801.   6.57636e-14,
  802.  -1.81807e-14,
  803.   4.6884e-15
  804. };
  805. static cheb_series fd_3half_d_cs = {
  806.   fd_3half_d_data,
  807.   24,
  808.   -1, 1,
  809.   16
  810. };
  811.  
  812.  
  813. /* Goano's modification of the Levin-u implementation.
  814.  * This is a simplification of the original WHIZ algorithm.
  815.  * See [Fessler et al., ACM Toms 9, 346 (1983)].
  816.  */
  817. static
  818. int
  819. fd_whiz(const double term, const int iterm,
  820.         double * qnum, double * qden,
  821.         double * result, double * s)
  822. {
  823.   if(iterm == 0) *s = 0.0;
  824.  
  825.   *s += term;
  826.  
  827.   qden[iterm] = 1.0/(term*(iterm+1.0)*(iterm+1.0));
  828.   qnum[iterm] = *s * qden[iterm];
  829.  
  830.   if(iterm > 0) {
  831.     double factor = 1.0;
  832.     double ratio  = iterm/(iterm+1.0);
  833.     int j;
  834.     for(j=iterm-1; j>=0; j--) {
  835.       double c = factor * (j+1.0) / (iterm+1.0);
  836.       factor *= ratio;
  837.       qden[j] = qden[j+1] - c * qden[j];
  838.       qnum[j] = qnum[j+1] - c * qnum[j];
  839.     }
  840.   }
  841.  
  842.   *result = qnum[0] / qden[0];
  843.   return GSL_SUCCESS;
  844. }
  845.  
  846.  
  847. /* Handle case of integer j <= -2.
  848.  */
  849. static
  850. int
  851. fd_nint(const int j, const double x, gsl_sf_result * result)
  852. {
  853. /*    const int nsize = 100 + 1; */
  854.   enum {
  855.     nsize = 100+1
  856.   };
  857.   double qcoeff[nsize];
  858.  
  859.   if(j >= -1) {
  860.     result->val = 0.0;
  861.     result->err = 0.0;
  862.     GSL_ERROR ("error", GSL_ESANITY);
  863.   }
  864.   else if(j < -(nsize)) {
  865.     result->val = 0.0;
  866.     result->err = 0.0;
  867.     GSL_ERROR ("error", GSL_EUNIMPL);
  868.   }
  869.   else {
  870.     double a, p, f;
  871.     int i, k;
  872.     int n = -(j+1);
  873.  
  874.     qcoeff[1] = 1.0;
  875.  
  876.     for(k=2; k<=n; k++) {
  877.       qcoeff[k] = -qcoeff[k-1];
  878.       for(i=k-1; i>=2; i--) {
  879.         qcoeff[i] = i*qcoeff[i] - (k-(i-1))*qcoeff[i-1];
  880.       }
  881.     }
  882.  
  883.     if(x >= 0.0) {
  884.       a = exp(-x);
  885.       f = qcoeff[1];
  886.       for(i=2; i<=n; i++) {
  887.         f = f*a + qcoeff[i];
  888.       }
  889.     }
  890.     else {
  891.       a = exp(x);
  892.       f = qcoeff[n];
  893.       for(i=n-1; i>=1; i--) {
  894.         f = f*a + qcoeff[i];
  895.       }
  896.     }
  897.  
  898.     p = gsl_sf_pow_int(1.0+a, j);
  899.     result->val = f*a*p;
  900.     result->err = 3.0 * GSL_DBL_EPSILON * fabs(f*a*p);
  901.     return GSL_SUCCESS;
  902.   }
  903. }
  904.  
  905.  
  906. /* x < 0
  907.  */
  908. static
  909. int
  910. fd_neg(const double j, const double x, gsl_sf_result * result)
  911. {
  912.   enum {
  913.     itmax = 100,
  914.     qsize = 100+1
  915.   };
  916. /*    const int itmax = 100; */
  917. /*    const int qsize = 100 + 1; */
  918.   double qnum[qsize], qden[qsize];
  919.  
  920.   if(x < GSL_LOG_DBL_MIN) {
  921.     result->val = 0.0;
  922.     result->err = 0.0;
  923.     return GSL_SUCCESS;
  924.   }
  925.   else if(x < -1.0 && x < -fabs(j+1.0)) {
  926.     /* Simple series implementation. Avoid the
  927.      * complexity and extra work of the series
  928.      * acceleration method below.
  929.      */
  930.     double ex   = exp(x);
  931.     double term = ex;
  932.     double sum  = term;
  933.     int n;
  934.     for(n=2; n<100; n++) {
  935.       double rat = (n-1.0)/n;
  936.       double p   = pow(rat, j+1.0);
  937.       term *= -ex * p;
  938.       sum  += term;
  939.       if(fabs(term/sum) < GSL_DBL_EPSILON) break;
  940.     }
  941.     result->val = sum;
  942.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(sum);
  943.     return GSL_SUCCESS;
  944.   }
  945.   else {
  946.     double s;
  947.     double xn = x;
  948.     double ex  = -exp(x);
  949.     double enx = -ex;
  950.     double f = 0.0;
  951.     double f_previous;
  952.     int jterm;
  953.     for(jterm=0; jterm<=itmax; jterm++) {
  954.       double p = pow(jterm+1.0, j+1.0);
  955.       double term = enx/p;
  956.       f_previous = f;
  957.       fd_whiz(term, jterm, qnum, qden, &f, &s);
  958.       xn += x;
  959.       if(fabs(f-f_previous) < fabs(f)*2.0*GSL_DBL_EPSILON || xn < GSL_LOG_DBL_MIN) break;
  960.       enx *= ex;
  961.     }
  962.  
  963.     result->val  = f;
  964.     result->err  = fabs(f-f_previous);
  965.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(f);
  966.  
  967.     if(jterm == itmax)
  968.       GSL_ERROR ("error", GSL_EMAXITER);
  969.     else
  970.       return GSL_SUCCESS;
  971.   }
  972. }
  973.  
  974.  
  975. /* asymptotic expansion
  976.  * j + 2.0 > 0.0
  977.  */
  978. static
  979. int
  980. fd_asymp(const double j, const double x, gsl_sf_result * result)
  981. {
  982.   const int j_integer = ( fabs(j - floor(j+0.5)) < 100.0*GSL_DBL_EPSILON );
  983.   const int itmax = 200;
  984.   gsl_sf_result lg;
  985.   int stat_lg = gsl_sf_lngamma_e(j + 2.0, &lg);
  986.   double seqn_val = 0.5;
  987.   double seqn_err = 0.0;
  988.   double xm2  = (1.0/x)/x;
  989.   double xgam = 1.0;
  990.   double add = GSL_DBL_MAX;
  991.   double cos_term;
  992.   double ln_x;
  993.   double ex_term_1;
  994.   double ex_term_2;
  995.   gsl_sf_result fneg;
  996.   gsl_sf_result ex_arg;
  997.   gsl_sf_result ex;
  998.   int stat_fneg;
  999.   int stat_e;
  1000.   int n;
  1001.   for(n=1; n<=itmax; n++) {
  1002.     double add_previous = add;
  1003.     gsl_sf_result eta;
  1004.     gsl_sf_eta_int_e(2*n, &eta);
  1005.     xgam = xgam * xm2 * (j + 1.0 - (2*n-2)) * (j + 1.0 - (2*n-1));
  1006.     add  = eta.val * xgam;
  1007.     if(!j_integer && fabs(add) > fabs(add_previous)) break;
  1008.     if(fabs(add/seqn_val) < GSL_DBL_EPSILON) break;
  1009.     seqn_val += add;
  1010.     seqn_err += 2.0 * GSL_DBL_EPSILON * fabs(add);
  1011.   }
  1012.   seqn_err += fabs(add);
  1013.  
  1014.   stat_fneg = fd_neg(j, -x, &fneg);
  1015.   ln_x = log(x);
  1016.   ex_term_1 = (j+1.0)*ln_x;
  1017.   ex_term_2 = lg.val;
  1018.   ex_arg.val = ex_term_1 - ex_term_2; /*(j+1.0)*ln_x - lg.val; */
  1019.   ex_arg.err = GSL_DBL_EPSILON*(fabs(ex_term_1) + fabs(ex_term_2)) + lg.err;
  1020.   stat_e    = gsl_sf_exp_err_e(ex_arg.val, ex_arg.err, &ex);
  1021.   cos_term  = cos(j*M_PI);
  1022.   result->val  = cos_term * fneg.val + 2.0 * seqn_val * ex.val;
  1023.   result->err  = fabs(2.0 * ex.err * seqn_val);
  1024.   result->err += fabs(2.0 * ex.val * seqn_err);
  1025.   result->err += fabs(cos_term) * fneg.err;
  1026.   result->err += 4.0 * GSL_DBL_EPSILON * fabs(result->val);
  1027.   return GSL_ERROR_SELECT_3(stat_e, stat_fneg, stat_lg);
  1028. }
  1029.  
  1030.  
  1031. /* Series evaluation for small x, generic j.
  1032.  * [Goano (8)]
  1033.  */
  1034. #if 0
  1035. static
  1036. int
  1037. fd_series(const double j, const double x, double * result)
  1038. {
  1039.   const int nmax = 1000;
  1040.   int n;
  1041.   double sum = 0.0;
  1042.   double prev;
  1043.   double pow_factor = 1.0;
  1044.   double eta_factor;
  1045.   gsl_sf_eta_e(j + 1.0, &eta_factor);
  1046.   prev = pow_factor * eta_factor;
  1047.   sum += prev;
  1048.   for(n=1; n<nmax; n++) {
  1049.     double term;
  1050.     gsl_sf_eta_e(j+1.0-n, &eta_factor);
  1051.     pow_factor *= x/n;
  1052.     term = pow_factor * eta_factor;
  1053.     sum += term;
  1054.     if(fabs(term/sum) < GSL_DBL_EPSILON && fabs(prev/sum) < GSL_DBL_EPSILON) break;
  1055.     prev = term;
  1056.   }
  1057.  
  1058.   *result = sum;
  1059.   return GSL_SUCCESS;
  1060. }
  1061. #endif /* 0 */
  1062.  
  1063.  
  1064. /* Series evaluation for small x > 0, integer j > 0; x < Pi.
  1065.  * [Goano (8)]
  1066.  */
  1067. static
  1068. int
  1069. fd_series_int(const int j, const double x, gsl_sf_result * result)
  1070. {
  1071.   int n;
  1072.   double sum = 0.0;
  1073.   double del;
  1074.   double pow_factor = 1.0;
  1075.   gsl_sf_result eta_factor;
  1076.   gsl_sf_eta_int_e(j + 1, &eta_factor);
  1077.   del  = pow_factor * eta_factor.val;
  1078.   sum += del;
  1079.  
  1080.   /* Sum terms where the argument
  1081.    * of eta() is positive.
  1082.    */
  1083.   for(n=1; n<=j+2; n++) {
  1084.     gsl_sf_eta_int_e(j+1-n, &eta_factor);
  1085.     pow_factor *= x/n;
  1086.     del  = pow_factor * eta_factor.val;
  1087.     sum += del;
  1088.     if(fabs(del/sum) < GSL_DBL_EPSILON) break;
  1089.   }
  1090.  
  1091.   /* Now sum the terms where eta() is negative.
  1092.    * The argument of eta() must be odd as well,
  1093.    * so it is convenient to transform the series
  1094.    * as follows:
  1095.    *
  1096.    * Sum[ eta(j+1-n) x^n / n!, {n,j+4,Infinity}]
  1097.    * = x^j / j! Sum[ eta(1-2m) x^(2m) j! / (2m+j)! , {m,2,Infinity}]
  1098.    *
  1099.    * We do not need to do this sum if j is large enough.
  1100.    */
  1101.   if(j < 32) {
  1102.     int m;
  1103.     gsl_sf_result jfact;
  1104.     double sum2;
  1105.     double pre2;
  1106.  
  1107.     gsl_sf_fact_e((unsigned int)j, &jfact);
  1108.     pre2 = gsl_sf_pow_int(x, j) / jfact.val;
  1109.  
  1110.     gsl_sf_eta_int_e(-3, &eta_factor);
  1111.     pow_factor = x*x*x*x / ((j+4)*(j+3)*(j+2)*(j+1));
  1112.     sum2 = eta_factor.val * pow_factor;
  1113.  
  1114.     for(m=3; m<24; m++) {
  1115.       gsl_sf_eta_int_e(1-2*m, &eta_factor);
  1116.       pow_factor *= x*x / ((j+2*m)*(j+2*m-1));
  1117.       sum2 += eta_factor.val * pow_factor;
  1118.     }
  1119.  
  1120.     sum += pre2 * sum2;
  1121.   }
  1122.  
  1123.   result->val = sum;
  1124.   result->err = 2.0 * GSL_DBL_EPSILON * fabs(sum);
  1125.  
  1126.   return GSL_SUCCESS;
  1127. }
  1128.  
  1129.  
  1130. /* series of hypergeometric functions for integer j > 0, x > 0
  1131.  * [Goano (7)]
  1132.  */
  1133. static
  1134. int
  1135. fd_UMseries_int(const int j, const double x, gsl_sf_result * result)
  1136. {
  1137.   const int nmax = 2000;
  1138.   double pre;
  1139.   double lnpre_val;
  1140.   double lnpre_err;
  1141.   double sum_even_val = 1.0;
  1142.   double sum_even_err = 0.0;
  1143.   double sum_odd_val  = 0.0;
  1144.   double sum_odd_err  = 0.0;
  1145.   int stat_sum;
  1146.   int stat_e;
  1147.   int stat_h = GSL_SUCCESS;
  1148.   int n;
  1149.  
  1150.   if(x < 500.0 && j < 80) {
  1151.     double p = gsl_sf_pow_int(x, j+1);
  1152.     gsl_sf_result g;
  1153.     gsl_sf_fact_e(j+1, &g); /* Gamma(j+2) */
  1154.     lnpre_val = 0.0;
  1155.     lnpre_err = 0.0;
  1156.     pre   = p/g.val;
  1157.   }
  1158.   else {
  1159.     double lnx = log(x);
  1160.     gsl_sf_result lg;
  1161.     gsl_sf_lngamma_e(j + 2.0, &lg);
  1162.     lnpre_val = (j+1.0)*lnx - lg.val;
  1163.     lnpre_err = 2.0 * GSL_DBL_EPSILON * fabs((j+1.0)*lnx) + lg.err;
  1164.     pre = 1.0;
  1165.   }
  1166.  
  1167.   /* Add up the odd terms of the sum.
  1168.    */
  1169.   for(n=1; n<nmax; n+=2) {
  1170.     double del_val;
  1171.     double del_err;
  1172.     gsl_sf_result U;
  1173.     gsl_sf_result M;
  1174.     int stat_h_U = gsl_sf_hyperg_U_int_e(1, j+2, n*x, &U);
  1175.     int stat_h_F = gsl_sf_hyperg_1F1_int_e(1, j+2, -n*x, &M);
  1176.     stat_h = GSL_ERROR_SELECT_3(stat_h, stat_h_U, stat_h_F);
  1177.     del_val = ((j+1.0)*U.val - M.val);
  1178.     del_err = (fabs(j+1.0)*U.err + M.err);
  1179.     sum_odd_val += del_val;
  1180.     sum_odd_err += del_err;
  1181.     if(fabs(del_val/sum_odd_val) < GSL_DBL_EPSILON) break;
  1182.   }
  1183.  
  1184.   /* Add up the even terms of the sum.
  1185.    */
  1186.   for(n=2; n<nmax; n+=2) {
  1187.     double del_val;
  1188.     double del_err;
  1189.     gsl_sf_result U;
  1190.     gsl_sf_result M;
  1191.     int stat_h_U = gsl_sf_hyperg_U_int_e(1, j+2, n*x, &U);
  1192.     int stat_h_F = gsl_sf_hyperg_1F1_int_e(1, j+2, -n*x, &M);
  1193.     stat_h = GSL_ERROR_SELECT_3(stat_h, stat_h_U, stat_h_F);
  1194.     del_val = ((j+1.0)*U.val - M.val);
  1195.     del_err = (fabs(j+1.0)*U.err + M.err);
  1196.     sum_even_val -= del_val;
  1197.     sum_even_err += del_err;
  1198.     if(fabs(del_val/sum_even_val) < GSL_DBL_EPSILON) break;
  1199.   }
  1200.  
  1201.   stat_sum = ( n >= nmax ? GSL_EMAXITER : GSL_SUCCESS );
  1202.   stat_e   = gsl_sf_exp_mult_err_e(lnpre_val, lnpre_err,
  1203.                                       pre*(sum_even_val + sum_odd_val),
  1204.                       pre*(sum_even_err + sum_odd_err),
  1205.                       result);
  1206.   result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  1207.  
  1208.   return GSL_ERROR_SELECT_3(stat_e, stat_h, stat_sum);
  1209. }
  1210.  
  1211.  
  1212. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  1213.  
  1214. /* [Goano (4)] */
  1215. int gsl_sf_fermi_dirac_m1_e(const double x, gsl_sf_result * result)
  1216. {
  1217.   if(x < GSL_LOG_DBL_MIN) {
  1218.     UNDERFLOW_ERROR(result);
  1219.   }
  1220.   else if(x < 0.0) {
  1221.     const double ex = exp(x);
  1222.     result->val = ex/(1.0+ex);
  1223.     result->err = 2.0 * (fabs(x) + 1.0) * GSL_DBL_EPSILON * fabs(result->val);
  1224.     return GSL_SUCCESS;
  1225.   }
  1226.   else {
  1227.     double ex = exp(-x);
  1228.     result->val = 1.0/(1.0 + ex);
  1229.     result->err = 2.0 * GSL_DBL_EPSILON * (x + 1.0) * ex;
  1230.     return GSL_SUCCESS;
  1231.   }
  1232. }
  1233.  
  1234.  
  1235. /* [Goano (3)] */
  1236. int gsl_sf_fermi_dirac_0_e(const double x, gsl_sf_result * result)
  1237. {
  1238.   if(x < GSL_LOG_DBL_MIN) {
  1239.     UNDERFLOW_ERROR(result);
  1240.   }
  1241.   else if(x < -5.0) {
  1242.     double ex  = exp(x);
  1243.     double ser = 1.0 - ex*(0.5 - ex*(1.0/3.0 - ex*(1.0/4.0 - ex*(1.0/5.0 - ex/6.0))));
  1244.     result->val = ex * ser;
  1245.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  1246.     return GSL_SUCCESS;
  1247.   }
  1248.   else if(x < 10.0) {
  1249.     result->val = log(1.0 + exp(x));
  1250.     result->err = fabs(x * GSL_DBL_EPSILON);
  1251.     return GSL_SUCCESS;
  1252.   }
  1253.   else {
  1254.     double ex = exp(-x);
  1255.     result->val = x + ex * (1.0 - 0.5*ex + ex*ex/3.0 - ex*ex*ex/4.0);
  1256.     result->err = (x + ex) * GSL_DBL_EPSILON;
  1257.     return GSL_SUCCESS;
  1258.   }
  1259. }
  1260.  
  1261.  
  1262. int gsl_sf_fermi_dirac_1_e(const double x, gsl_sf_result * result)
  1263. {
  1264.   if(x < GSL_LOG_DBL_MIN) {
  1265.     UNDERFLOW_ERROR(result);
  1266.   }
  1267.   else if(x < -1.0) {
  1268.     /* series [Goano (6)]
  1269.      */
  1270.     double ex   = exp(x);
  1271.     double term = ex;
  1272.     double sum  = term;
  1273.     int n;
  1274.     for(n=2; n<100 ; n++) {
  1275.       double rat = (n-1.0)/n;
  1276.       term *= -ex * rat * rat;
  1277.       sum  += term;
  1278.       if(fabs(term/sum) < GSL_DBL_EPSILON) break;
  1279.     }
  1280.     result->val = sum;
  1281.     result->err = 2.0 * fabs(sum) * GSL_DBL_EPSILON;
  1282.     return GSL_SUCCESS;
  1283.   }
  1284.   else if(x < 1.0) {
  1285.     return cheb_eval_e(&fd_1_a_cs, x, result);
  1286.   }
  1287.   else if(x < 4.0) {
  1288.     double t = 2.0/3.0*(x-1.0) - 1.0;
  1289.     return cheb_eval_e(&fd_1_b_cs, t, result);
  1290.   }
  1291.   else if(x < 10.0) {
  1292.     double t = 1.0/3.0*(x-4.0) - 1.0;
  1293.     return cheb_eval_e(&fd_1_c_cs, t, result);
  1294.   }
  1295.   else if(x < 30.0) {
  1296.     double t = 0.1*x - 2.0;
  1297.     gsl_sf_result c;
  1298.     cheb_eval_e(&fd_1_d_cs, t, &c);
  1299.     result->val  = c.val * x*x;
  1300.     result->err  = c.err * x*x + GSL_DBL_EPSILON * fabs(result->val);
  1301.     return GSL_SUCCESS;
  1302.   }
  1303.   else if(x < 1.0/GSL_SQRT_DBL_EPSILON) {
  1304.     double t = 60.0/x - 1.0;
  1305.     gsl_sf_result c;
  1306.     cheb_eval_e(&fd_1_e_cs, t, &c);
  1307.     result->val  = c.val * x*x;
  1308.     result->err  = c.err * x*x + GSL_DBL_EPSILON * fabs(result->val);
  1309.     return GSL_SUCCESS;
  1310.   }
  1311.   else if(x < GSL_SQRT_DBL_MAX) {
  1312.     result->val = 0.5 * x*x;
  1313.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  1314.     return GSL_SUCCESS;
  1315.   }
  1316.   else {
  1317.     OVERFLOW_ERROR(result);
  1318.   }
  1319. }
  1320.  
  1321.  
  1322. int gsl_sf_fermi_dirac_2_e(const double x, gsl_sf_result * result)
  1323. {
  1324.   if(x < GSL_LOG_DBL_MIN) {
  1325.     UNDERFLOW_ERROR(result);
  1326.   }
  1327.   else if(x < -1.0) {
  1328.     /* series [Goano (6)]
  1329.      */
  1330.     double ex   = exp(x);
  1331.     double term = ex;
  1332.     double sum  = term;
  1333.     int n;
  1334.     for(n=2; n<100 ; n++) {
  1335.       double rat = (n-1.0)/n;
  1336.       term *= -ex * rat * rat * rat;
  1337.       sum  += term;
  1338.       if(fabs(term/sum) < GSL_DBL_EPSILON) break;
  1339.     }
  1340.     result->val = sum;
  1341.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(sum);
  1342.     return GSL_SUCCESS;
  1343.   }
  1344.   else if(x < 1.0) {
  1345.     return cheb_eval_e(&fd_2_a_cs, x, result);
  1346.   }
  1347.   else if(x < 4.0) {
  1348.     double t = 2.0/3.0*(x-1.0) - 1.0;
  1349.     return cheb_eval_e(&fd_2_b_cs, t, result);
  1350.   }
  1351.   else if(x < 10.0) {
  1352.     double t = 1.0/3.0*(x-4.0) - 1.0;
  1353.     return cheb_eval_e(&fd_2_c_cs, t, result);
  1354.   }
  1355.   else if(x < 30.0) {
  1356.     double t = 0.1*x - 2.0;
  1357.     gsl_sf_result c;
  1358.     cheb_eval_e(&fd_2_d_cs, t, &c);
  1359.     result->val  = c.val * x*x*x;
  1360.     result->err  = c.err * x*x*x + 3.0 * GSL_DBL_EPSILON * fabs(result->val);
  1361.     return GSL_SUCCESS;
  1362.   }
  1363.   else if(x < 1.0/GSL_ROOT3_DBL_EPSILON) {
  1364.     double t = 60.0/x - 1.0;
  1365.     gsl_sf_result c;
  1366.     cheb_eval_e(&fd_2_e_cs, t, &c);
  1367.     result->val = c.val * x*x*x;
  1368.     result->err = c.err * x*x*x + 3.0 * GSL_DBL_EPSILON * fabs(result->val);
  1369.     return GSL_SUCCESS;
  1370.   }
  1371.   else if(x < GSL_ROOT3_DBL_MAX) {
  1372.     result->val = 1.0/6.0 * x*x*x;
  1373.     result->err = 3.0 * GSL_DBL_EPSILON * fabs(result->val);
  1374.     return GSL_SUCCESS;
  1375.   }
  1376.   else {
  1377.     OVERFLOW_ERROR(result);
  1378.   }
  1379. }
  1380.  
  1381.  
  1382. int gsl_sf_fermi_dirac_int_e(const int j, const double x, gsl_sf_result * result)
  1383. {
  1384.   if(j < -1) {
  1385.     return fd_nint(j, x, result);
  1386.   }
  1387.   else if (j == -1) {
  1388.     return gsl_sf_fermi_dirac_m1_e(x, result);
  1389.   }
  1390.   else if(j == 0) {
  1391.     return gsl_sf_fermi_dirac_0_e(x, result);
  1392.   }
  1393.   else if(j == 1) {
  1394.     return gsl_sf_fermi_dirac_1_e(x, result);
  1395.   }
  1396.   else if(j == 2) {
  1397.     return gsl_sf_fermi_dirac_2_e(x, result);
  1398.   }
  1399.   else if(x < 0.0) {
  1400.     return fd_neg(j, x, result);
  1401.   }
  1402.   else if(x == 0.0) {
  1403.     return gsl_sf_eta_int_e(j+1, result);
  1404.   }
  1405.   else if(x < 1.5) {
  1406.     return fd_series_int(j, x, result);
  1407.   }
  1408.   else {
  1409.     gsl_sf_result fasymp;
  1410.     int stat_asymp = fd_asymp(j, x, &fasymp);
  1411.  
  1412.     if(stat_asymp == GSL_SUCCESS) {
  1413.       result->val  = fasymp.val;
  1414.       result->err  = fasymp.err;
  1415.       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  1416.       return stat_asymp;
  1417.     }
  1418.     else {
  1419.       return fd_UMseries_int(j, x, result);
  1420.     }
  1421.   }
  1422. }
  1423.  
  1424.  
  1425. int gsl_sf_fermi_dirac_mhalf_e(const double x, gsl_sf_result * result)
  1426. {
  1427.   if(x < GSL_LOG_DBL_MIN) {
  1428.     UNDERFLOW_ERROR(result);
  1429.   }
  1430.   else if(x < -1.0) {
  1431.     /* series [Goano (6)]
  1432.      */
  1433.     double ex   = exp(x);
  1434.     double term = ex;
  1435.     double sum  = term;
  1436.     int n;
  1437.     for(n=2; n<200 ; n++) {
  1438.       double rat = (n-1.0)/n;
  1439.       term *= -ex * sqrt(rat);
  1440.       sum  += term;
  1441.       if(fabs(term/sum) < GSL_DBL_EPSILON) break;
  1442.     }
  1443.     result->val = sum;
  1444.     result->err = 2.0 * fabs(sum) * GSL_DBL_EPSILON;
  1445.     return GSL_SUCCESS;
  1446.   }
  1447.   else if(x < 1.0) {
  1448.     return cheb_eval_e(&fd_mhalf_a_cs, x, result);
  1449.   }
  1450.   else if(x < 4.0) {
  1451.     double t = 2.0/3.0*(x-1.0) - 1.0;
  1452.     return cheb_eval_e(&fd_mhalf_b_cs, t, result);
  1453.   }
  1454.   else if(x < 10.0) {
  1455.     double t = 1.0/3.0*(x-4.0) - 1.0;
  1456.     return cheb_eval_e(&fd_mhalf_c_cs, t, result);
  1457.   }
  1458.   else if(x < 30.0) {
  1459.     double rtx = sqrt(x);
  1460.     double t = 0.1*x - 2.0;
  1461.     gsl_sf_result c;
  1462.     cheb_eval_e(&fd_mhalf_d_cs, t, &c);
  1463.     result->val  = c.val * rtx;
  1464.     result->err  = c.err * rtx + 0.5 * GSL_DBL_EPSILON * fabs(result->val);
  1465.     return GSL_SUCCESS;
  1466.   }
  1467.   else {
  1468.     return fd_asymp(-0.5, x, result);
  1469.   }
  1470. }
  1471.  
  1472.  
  1473. int gsl_sf_fermi_dirac_half_e(const double x, gsl_sf_result * result)
  1474. {
  1475.   if(x < GSL_LOG_DBL_MIN) {
  1476.     UNDERFLOW_ERROR(result);
  1477.   }
  1478.   else if(x < -1.0) {
  1479.     /* series [Goano (6)]
  1480.      */
  1481.     double ex   = exp(x);
  1482.     double term = ex;
  1483.     double sum  = term;
  1484.     int n;
  1485.     for(n=2; n<100 ; n++) {
  1486.       double rat = (n-1.0)/n;
  1487.       term *= -ex * rat * sqrt(rat);
  1488.       sum  += term;
  1489.       if(fabs(term/sum) < GSL_DBL_EPSILON) break;
  1490.     }
  1491.     result->val = sum;
  1492.     result->err = 2.0 * fabs(sum) * GSL_DBL_EPSILON;
  1493.     return GSL_SUCCESS;
  1494.   }
  1495.   else if(x < 1.0) {
  1496.     return cheb_eval_e(&fd_half_a_cs, x, result);
  1497.   }
  1498.   else if(x < 4.0) {
  1499.     double t = 2.0/3.0*(x-1.0) - 1.0;
  1500.     return cheb_eval_e(&fd_half_b_cs, t, result);
  1501.   }
  1502.   else if(x < 10.0) {
  1503.     double t = 1.0/3.0*(x-4.0) - 1.0;
  1504.     return cheb_eval_e(&fd_half_c_cs, t, result);
  1505.   }
  1506.   else if(x < 30.0) {
  1507.     double x32 = x*sqrt(x);
  1508.     double t = 0.1*x - 2.0;
  1509.     gsl_sf_result c;
  1510.     cheb_eval_e(&fd_half_d_cs, t, &c);
  1511.     result->val = c.val * x32;
  1512.     result->err = c.err * x32 + 1.5 * GSL_DBL_EPSILON * fabs(result->val);
  1513.     return GSL_SUCCESS;
  1514.   }
  1515.   else {
  1516.     return fd_asymp(0.5, x, result);
  1517.   }
  1518. }
  1519.  
  1520.  
  1521. int gsl_sf_fermi_dirac_3half_e(const double x, gsl_sf_result * result)
  1522. {
  1523.   if(x < GSL_LOG_DBL_MIN) {
  1524.     UNDERFLOW_ERROR(result);
  1525.   }
  1526.   else if(x < -1.0) {
  1527.     /* series [Goano (6)]
  1528.      */
  1529.     double ex   = exp(x);
  1530.     double term = ex;
  1531.     double sum  = term;
  1532.     int n;
  1533.     for(n=2; n<100 ; n++) {
  1534.       double rat = (n-1.0)/n;
  1535.       term *= -ex * rat * rat * sqrt(rat);
  1536.       sum  += term;
  1537.       if(fabs(term/sum) < GSL_DBL_EPSILON) break;
  1538.     }
  1539.     result->val = sum;
  1540.     result->err = 2.0 * fabs(sum) * GSL_DBL_EPSILON;
  1541.     return GSL_SUCCESS;
  1542.   }
  1543.   else if(x < 1.0) {
  1544.     return cheb_eval_e(&fd_3half_a_cs, x, result);
  1545.   }
  1546.   else if(x < 4.0) {
  1547.     double t = 2.0/3.0*(x-1.0) - 1.0;
  1548.     return cheb_eval_e(&fd_3half_b_cs, t, result);
  1549.   }
  1550.   else if(x < 10.0) {
  1551.     double t = 1.0/3.0*(x-4.0) - 1.0;
  1552.     return cheb_eval_e(&fd_3half_c_cs, t, result);
  1553.   }
  1554.   else if(x < 30.0) {
  1555.     double x52 = x*x*sqrt(x);
  1556.     double t = 0.1*x - 2.0;
  1557.     gsl_sf_result c;
  1558.     cheb_eval_e(&fd_3half_d_cs, t, &c);
  1559.     result->val = c.val * x52;
  1560.     result->err = c.err * x52 + 2.5 * GSL_DBL_EPSILON * fabs(result->val);
  1561.     return GSL_SUCCESS;
  1562.   }
  1563.   else {
  1564.     return fd_asymp(1.5, x, result);
  1565.   }
  1566. }
  1567.  
  1568. /* [Goano p. 222] */
  1569. int gsl_sf_fermi_dirac_inc_0_e(const double x, const double b, gsl_sf_result * result)
  1570. {
  1571.   if(b < 0.0) {
  1572.     DOMAIN_ERROR(result);
  1573.   }
  1574.   else {
  1575.     double arg = b - x;
  1576.     gsl_sf_result f0;
  1577.     int status = gsl_sf_fermi_dirac_0_e(arg, &f0);
  1578.     result->val = f0.val - arg;
  1579.     result->err = f0.err + GSL_DBL_EPSILON * (fabs(x) + fabs(b));
  1580.     return status;
  1581.   }
  1582. }
  1583.  
  1584.  
  1585.  
  1586. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  1587.  
  1588. #include "eval.h"
  1589.  
  1590. double gsl_sf_fermi_dirac_m1(const double x)
  1591. {
  1592.   EVAL_RESULT(gsl_sf_fermi_dirac_m1_e(x, &result));
  1593. }
  1594.  
  1595. double gsl_sf_fermi_dirac_0(const double x)
  1596. {
  1597.   EVAL_RESULT(gsl_sf_fermi_dirac_0_e(x, &result));
  1598. }
  1599.  
  1600. double gsl_sf_fermi_dirac_1(const double x)
  1601. {
  1602.   EVAL_RESULT(gsl_sf_fermi_dirac_1_e(x, &result));
  1603. }
  1604.  
  1605. double gsl_sf_fermi_dirac_2(const double x)
  1606. {
  1607.   EVAL_RESULT(gsl_sf_fermi_dirac_2_e(x, &result));
  1608. }
  1609.  
  1610. double gsl_sf_fermi_dirac_int(const int j, const double x)
  1611. {
  1612.   EVAL_RESULT(gsl_sf_fermi_dirac_int_e(j, x, &result));
  1613. }
  1614.  
  1615. double gsl_sf_fermi_dirac_mhalf(const double x)
  1616. {
  1617.   EVAL_RESULT(gsl_sf_fermi_dirac_mhalf_e(x, &result));
  1618. }
  1619.  
  1620. double gsl_sf_fermi_dirac_half(const double x)
  1621. {
  1622.   EVAL_RESULT(gsl_sf_fermi_dirac_half_e(x, &result));
  1623. }
  1624.  
  1625. double gsl_sf_fermi_dirac_3half(const double x)
  1626. {
  1627.   EVAL_RESULT(gsl_sf_fermi_dirac_3half_e(x, &result));
  1628. }
  1629.  
  1630. double gsl_sf_fermi_dirac_inc_0(const double x, const double b)
  1631. {
  1632.   EVAL_RESULT(gsl_sf_fermi_dirac_inc_0_e(x, b, &result));
  1633. }
  1634.